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ABSTRACT 

Delaunay- Voronoi  mesh  systems  provide  a  generalization  of  the  classical  rectangular  stag¬ 
gered  meshes  to  unstructured  meshes.  In  this  work,  it  is  shown  how  such  “covolume”  dis¬ 
cretizations  may  be  applied  to  div-curl  systems  in  three  dimensions.  Error  estimates  are 
proved  and  confirmed  by  a  numerical  illustration. 
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1  Introduction.  This  paper  deals  with  the  problem  of  numerically  determining 
a  three  dimensional  vector  field  from  its  divergence  and  curl.  The  need  to  do  this 
occurs  in  several  fields  including  fluid  dynamics  and  electromagnetics.  The  governing 
equations  are  often  solved  using  indirect  approaches  including  potential  formulations 
or  Biot-Savart  type  integrals  [5]  and  least  squares  [1]. 

A  good  reason  to  be  cautious  about  direct  discretizations  is  that  the  equations  are 
overdetermined,  being  a  system  of  four  equations  in  three  unknowns,  and  it  is  not 
immediately  evident  what  effect  this  will  have  at  the  discrete  level.  Nevertheless, 
since  potential  formulations  can  suffer  from  spurious  mode  problems  [6]  and  the 
Biot-Savart  approach  requires  special  handling  of  boundaries,  a  direct  treatment 
may  be  preferable.  We  present  such  a  treatment  in  this  paper. 

The  method  we  discuss  is  a  generalization  of  that  presented  in  [3]  for  planar  systems. 
As  in  the  planar  case  we  use  dual  mesh  systems  made  up  of  “complementary  vol¬ 
umes”  (or  “covolumes”)  of  which  the  prime  examples  are  classical  Voronoi-Delaunay 
mesh  pairs.  In  the  Voronoi-Delaunay  approach  the  domain  is  suitably  partitioned 
into  tetrahedra  which  are  considered  to  form  a  primal  (Delaunay)  mesh  in  the  do¬ 
main  D.  A  dual  (Voronoi)  mesh  is  formed  by  connecting  the  circumcenters  of  the 
tetrahedra,  giving  a  set  of  convex  polygons  each  of  which  encloses  a  unique  tetrahe¬ 
dral  vertex.  Combinatorial  and  topological  properties  of  the  primal  and  dual  mesh 
systems  play  an  important  part  in  our  analysis,  allowing  us  to  conclude,  for  example, 
that  the  solution  of  the  discretized  equations  is  exactly  determined  by  the  data  if  the 
same  is  true  of  the  original  system.  On  the  other  hand,  the  metric  properties  of  the 
mesh  system  are  used  in  deriving  the  error  estimates.  The  interplay  of  topological 
and  metric  properties  is  a  common  feature  of  covolume  discretizations. 

The  paper  contains  six  sections.  The  second  section  deals  with  combinatorial  prop¬ 
erties  of  the  mesh  system.  The  third  section  shows  how  discrete  vector  fields  are 
defined  on  the  meshes  and  defines  the  basic  integral  operations  for  the  discrete 
fields.  These  definitions  are  set  up  so  that  an  important  orthogonahty  property  of 
vector  fields  is  preserved  (theorem  4).  The  discrete  approximations  and  the  basic 
error  estimates  are  in  section  4.  In  section  -5  we  specialize  the  results  to  rectangu¬ 
lar  meshes  for  which  we  are  able  to  obtain  improved  convergence  rates.  The  final 
section  contains  some  basic  numerical  verifications. 


2  Mesh  notations  and  some  basic  results.  We  will  prove  our  results  for  the 
case  of  a  domain  D  which  is  a  bounded  polyhedron  in  with  boundary  T.  The 
results  we  present  are  generally  valid  for  multiply  connected  domains  and  domains 
with  cutouts.  The  techniques  for  extending  the  results  are  similar  to  those  used  in  [3] 
for  the  two  dimensional  case  although  there  is  greater  variety  in  three  dimensional 
cases  depending  on  the  genus  of  the  surface(s)  bounding  the  obstacle(s). 
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Assume  that  U  has  a  primal  family  of  finite  element  style  tetrahedral  partitions, 

parameterized  by  the  maximum  side  length  which  is  generlcaUy  denoted  by  h.  We 
wiU  assume  that  the  ratios  of  the  radii  of  the  circumscribing  spheres  and  the  in¬ 
scribed  spheres  of  all  the  individual  tetrahedra  are  uniformly  bounded  above  and 
below  as  h  approaches  0.  A  dual  mesh  is  formed  by  connecting  adjacent  tetrahedral 
circumcenters  and,  in  the  case  of  tetrahedra  with  a  face  on  a  boundary,  by  connect¬ 
ing  their  circumcenters  with  those  of  their  boundary  faces.  By  elementary  geometry 
the  connecting  edges  are  perpendicular  to  the  associated  tetrahedral  faces.  These 
connections  also  form  the  edges  of  a  set  of  polyhedra.  It  follows  from  elementary- 
geometry  that  the  edges  of  the  tetrahedra  are  perpendicular  to  and  in  one-to-one 
correspondence  with  the  faces  of  the  dual  polyhedra  (covolumes)  and  dually.  The 
reciprocal  orthogonality  between  edges  and  faces  is  the  key  to  the  results  which 
follow.  A  survey  of  algorithms  for  generating  Voronoi-Delaunay  meshes  is  given 
in  [4]. 

For  an  entirely  arbitrary  tetrahedral  partition  the  structure  of  the  covolume  partition 
can  be  complicated.  To  obtain  a  practical  algorithm  we  will  restrict  the  tetrahedral 
partition  so  that  the  dual  partition  has  two  (possibly  dependent)  properties.  These 
properties  are  (1)  interior  covolumes  have  faces  which  are  simple  (planar)  polygons 
and  (2)  each  tetrahedral  Interior  node  lies  in  exactly  one  covolume.  This  does  restrict 
the  primal  tesselation  somewhat  but  not  unduly.  For  instance  the  two  conditions 
are  certainly  satisfied  for  a  tetrahedral  partition  in  which  the  dihedral  angles  are  aU 
acute.  In  that  case,  the  two  meshes  form  a  Delaunay-Voronoi  pair  and  the  interior 
covolumes  are  convex.  In  addition  the  faces  of  the  covolumes  are  not  merely  simple 
but  convex  as  well.  More  generally  any  Delaunay-Voronoi  pair  will  satisfy  (1)  and 
(2). 

The  N  nodes  of  the  tetrahedral  mesh  are  assumed  to  be  numbered  sequentially 
in  some  convenient  way,  and  likewise  the  T  nodes  of  the  dual  mesh.  Similarly, 
the  F  faces  (edges)  and  E  edges  (faces)  of  the  primal  (dual)  mesh  are  sequentially 
numbered.  A  subscript  1  on  one  of  these  quantities  denotes  the  corresponding 
interior  number.  For  example  Fi  denotes  the  number  of  tetrahedral  interior  faces. 
The  individual  tetrahedra,  faces,  edges  and  nodes  of  the  primal  mesh  are  denoted 
by  Tj,  12 j,  Gk  and  vi  respectively.  Those  of  the  covolume  mesh  are  denoted  by  primed 
quantities  such  as  cr'.  A  direction  is  assigned  to  each  primal  edge  by  the  rule  that 
the  positive  direction  is  from  low  to  high  node  number.  The  dual  edges  are  directed 
by  the  corresponding  rule. 

We  will  use  the  following  standard  results: 

Lemma  1  Let  iV,  F,  E  and  T  denote  the  number  of  nodes,  faces,  edges  and  tetra¬ 
hedra  of  a  given  tesselation  of  a  polyhedral  domain.  Then 


N  -i-  F 

Ai  +  Fi 
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F  +  r  + 1 
Fi  +  r-  1. 


(1) 

(2) 


These  can  easily  be  proven  by  successive  deletions  of  tetrahedra  from  the  triangu¬ 
lation.  ■ 


3  Discrete  vector  fields. The  main  theorems  in  this  section  (theorems  3  and 
4)  are  discrete  analogs  of  theorems  in  vector  field  theory.  Theorem  3  asserts  the 
existence  of  a  scalar  potential  when  discrete  circulation  equations  are  satisfied  and 
theorem  4  gives  a  discrete  Helmholtz  decomposition  of  a  field  of  normal  components. 

Each  boundary  dfij  is  assumed  to  be  oriented  by  a  right  hand  rule  applied 

relative  to  the  directed  edge  cr'  (cr^). 

For  each  strictly  interior  dual  edge  a'-  we  can  form  a  vector  whose  kth  component  is 
the  sign  of  the  orientation  of  the  edge  relative  to  the  orientation  of  the  kth  strictly 
interior  dual  face.  From  these  vectors  we  obtain  the  Fi  x  E\  matrix  G  defined  as 
follows: 

{+1  if  cr'  is  oriented  positively  along 
-1  if  a'j  is  oriented  negatively  along  djji'f. 

0  if  a'-  does  not  meet  dfi').. 


G  has  rank  Ei  -  iV’i.  To  see  the  reason  for  this,  note  first  that  each  column  of  G  is 
associated  with  a  dual  face.  Now  take  any  interior  covolume  and  sum  the  associated 
columns,  each  weighted  ±1  according  to  whether  its  normal  points  out  of  or  into  the 
covolume.  Then  the  contributions  from  each  dual  edge  appear  twice  in  each  sum 
with  opposite  signs.  In  this  way  we  can  form  one  nuU  vector  of  G*  for  each  interior 
covolume.  It  is  clear  that  there  are  no  other  vectors  in  the  nuU  space  so  the  result 
follows. 


In  addition  to  G  we  will  use  another  matrix  B\  containing  the  orientations  of  the 
tetrahedral  faces  relative  to  the  outer  normals  of  the  tetrahedra.  This  matrix  is 
defined  for  interior  faces  and  has  dimensions  iq  x  T.  The  definition  of  Bi  is: 


(^1  )ji 


+  1  if  fij  is  oriented  with  dri 
<  —1  if  jj,j  is  oriented  against  dTi  , 

0  if  fij  does  not  meet  dvi. 


where  “oriented  with”  means  that  the  normal  direction  of  fij  is  parallel  to  the  outer 
normal  of  dr^.  The  T  dimensional  vector  with  I  in  every  position  is  in  the  nuU  space 
of  Bi  and  is  the  only  such  vector,  so  that  the  rank  of  H  is  T  —  1. 

In  addition  to  Bi  we  wiU  also  use  the  matrix  B  which  is  obtained  from  Bi  by  aug¬ 
menting  it  with  rows  corresponding  to  the  dual  edges  which  are  normal  to  boundary 
faces.  Recall  that  these  edges  pass  through  the  circumcenters  of  the  boundary  faces. 
The  additional  values  for  the  domain  of  B  are  associated  with  these  circumcenters. 


i?i  and  G  are  related  through  G^Bi  =  0.  To  see  this,  take  a  standard  basis  vector 
and  multiply  it  by  B.  Then  values  ±1  or  0  result  on  the  dual  edges  the  positive 
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(negative)  sign  being  associated  with  a  dual  edge  directed  out  of  (into)  the  node 

which  carries  the  unit  value.  Multiplication  by  G*  corresponds  to  forming  a  signed 
sum  of  these  values  around  the  dual  faces  where  the  signs  are  precisely  such  as  to 
effect  a  canceUation  of  the  nonzero  contributions.  Transposing  this  result  shows 
that  R{G)  C  N{B\)  where  N{-)  and  R[-)  denote  the  nuU  space  and  range  of  their 
arguments. 

Theorem  1  Let  v  G  R^'^  ■  Then  1  0  e  R^  such  that  v  -  Bi(f)  if  Gh  =  0. 

Proof.  Since  Bi  is  Fi  x  T  and  of  rank  T  -  1  and  using  lemma  1  we  have 

d\mN{B\ )  =  Fi-T+1 

=  Er-Ni. 


RecaUing  that  R(G)  C  N{B{)  and  since  dimi2(G)  =  £^i  -  iVi  we  have  N{B\)  = 
R(G).  Solvability  of  the  equation 


Bicj)  =  V 


holds  if 


(v.2)  =  0  'izeN{Bl), 


where  ( ,  )  denotes  the  standard  EucUdean  inner  product.  From  above  this  is  equiv¬ 
alent  to 


(tuGVO  =  0  \f^beR^'^ 


and  these  are  equivalent  to  the  theorem.  ■ 

The  covolume  algorithm  for  the  div-curl  equations  provides  approximations  to  the 
quantities  u-n,  where  u  denotes  the  exact  solution  vector  at  the  center  of  a  primal 
face  and  n  denotes  a  normal  to  the  face.  The  corresponding  approximate  quantities 
are  denoted  by  uj  where  j  indexes  the  faces. 

Any  set  of  normal  components  defined  on  faces  can  be  identified  with  R^ .  We  will 
introduce  an  inner  product  into  R^  by 

{u,v  fv  :=  Y. 

fij€Q 


where  Sj  is  the  area  of  fij  and  fi'  is  the  length  of  a'-.  The  resulting  inner  product 
space  is  denoted  by  ZY,  and  Uq  denotes  the  space 


ILq  :=  {u  eH:u  |r=  0}. 
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The  associated  norm  is  denoted  by  Clearly,  this  is  three  times  a  discrete 

norm. 

For  each  tetrahedron  r,  discrete  flux  and  divergence  operators  are  defined  on  U  by 

{bu),  := 

tijEdri 

and 

{Du)^  {bu)ilAi. 

By  Sj  we  mean  Sj  negatively  signed  if  the  corresponding  velocity  component  is 
directed  towards  the  inside  of  Tj  and  positively  signed  otherwise. 

For  each  interior  covolume  face  bk  discrete  circulation  and  curl  operators  are  defined 
by 

{Cu)k  :=  ufh'j 

and 

{Cu)k  :=  {Cu)k/s'i,. 

The  tilde  on  h'  means  that  this  quantity  is  to  be  taken  with  a  negative  sign  if  the 
dual  edge  is  directed  against  the  positive  sense  of  description  of  and  a  positive 
sign  otherwise. 

Denote  by  5  and  the  diagonal  matrices  S  :=  diag(sj)  and  :=  diag(h').  It 
can  be  checked  directly  that 

b  =  B^S 
C  =  G'iX', 

where  H'  denotes  the  restriction  of  Hi  to  interior  dual  edges.  We  also  define  differ¬ 
ence  operators 

Pb4> 

Pd>  :=  \/4>€R^. 

Corresponding  to  theorem  1  we  have 

Theorem  2  Let  v  G  R^'^  ■  Then  3  ^  G  R^  such  that  v  =  P(f)  if  Cv  =  0.1 

Now  introduce  two  subspaces  of  ILq: 

Zq  =  {u  G  Uq]  bu  =  0} 
kVo  =  {u  G  ZYo;  Cm  =  0}. 

Then  we  have 
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Theorem  3  Let  u  G  Zq,  v  G  Wq,  then 

(u,v)w  =  0. 

Proof.  From  theorem  2  and  using  summation  by  parts  we  have 

(u,  v)w  =  {u,  P4>)w  -  (w,  Pb(p)w  -  (Du,  (f>)  MS  €  . 

Since  u  G  2o  we  have  Du  =  0  and  the  theorem  follows.  ■ 


4  Div-curl  systems.  In  this  section  we  formulate  the  discrete  approximations 
for  div-curl  equations  and  show  that  they  have  a  unique  solution  (theorem  5).  The 
basic  error  estimate  is  given  in  theorem  6. 

We  consider  the  three  dimensional  div-curl  problem  for  u, 

divu  =  p  (3) 

curlu  =  uj  (4) 

u-n|r  =  /  .  (5) 

We  will  assume  that  the  following  compatibility  conditions  hold: 

divw  =  0  (6) 


/  pdv  =  [ 

Jq  Jr 


See  [2]  for  the  mathematical  background  to  this  problem. 
Equations  (3)-(5)  are  discretized  by 


where 


Up” 

i/  /* 

J  Hi 


Pj 


The  equations  (8)-(10)  are  not  all  independent.  In  fact,  there  are  at  least  Ei  relations 
between  the  equations  (9)  corresponding  to  the  fact  that  the  area  weighted  sums  of 
the  data  over  covolumes  is  zero  (corresponding  to  the  compatibility  condition  (7)) 
as  well  as  a  further  relation  which  is  equivalent  to  (7).  The  discrete  equations  do 
have  a  unique  solution  as  we  prove  next. 
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Theorem  4  Equaiions  (8)-(10)  have  a  unique  solution. 


Proof.  Consider  the  homogeneous  equations  with  =  0  and  /  =  0  in  (8)-(10). 
These  equations  have  at  least  one  solution  u,  and  this  solution  satisfies  u  e  Zq  and 
u  e  H-’o-  But  then  theorem  3  implies  m  =  0  which  implies  the  uniqueness  part  of  the 
result.  Let  K  denote  the  (T  +  iq  x  Fi  coefficient  matrix  of  the  linear  system;  then 
it  follows  that  K  has  rank  Fi- 

For  the  existence,  we  consider  the  linear  system  generated  by  (8)-(10).  and  let  K' 
denote  the  corresponding  augmented  matrix  of  K  of  order  (T  +  Ei)  x  (iq  +  1).  It 
follows  from  the  compatibility  conditions  that  this  matrix  has  rank  at  most  T  -  1  + 
El  —  N\  =  Fi  and  existence  follows  from  this.  ■ 

Next,  we  introduce  two  functions  and  u^^^  defined  by 

:=  —  /  u-ncl5  pj  G  fi 

V  I >  ^ 

If  (7j  is  connected  to  the  boundary  F,  u^^^  ;= 

Lemma  2  Assume  that  u  G  >  2,  then  we  have 

Proof.  Let  Hj  denote  the  polyhedron  associated  with  the  face  Hj  and  the  dual  edge 
cr'  which  is  obtained  by  connecting  the  circumcenters  of  the  primal  cells  which  share 
Pj  to  its  vertices.  By  a  Sobolev  embedding  theorem  we  have 

W^-p(kj)  ^  LHk])  p>2  r  =  l,2, 

where  k]  =  <t'-  and  nj  =  fij.  It  follows  that  V  :  u  ^  -  u^^^  defines  a  bounded 

linear  functional  on  W^’P{Kj).  Clearly,  is  zero  for  constant  u  so  we  have 


and  a  standard  rescaling  argument  then  shows  that  the  right  side  is 

where  C  is  independent  of  h  and  u.  Summing  over  the  domain  we  obtain 

\uf^  -  uf\h,h’^  <  c  5: 


<  Y.  HL..- 

IJ.J& 


( 


5  Uniform  meshes.  In  this  section  we  wiU  look  more  closely  at  the  situation  in 
which  the  domain  fi  is  a  cube  and  the  mesh  is  a  uniform  cubic  mesh.  For  this  case 
the  dual  mesh  is  also  uniform  cubic,  being  shifted  by  one  half  the  mesh  spacing  in 
each  coordinate  direction  from  the  primal  mesh. 

For  this  uniform  mesh  case  the  basic  error  estimate  can  be  improved.  The  difference 
in  the  uniform  mesh  case  arises  from  the  improved  approximation  theory  which  is 
possible.  Apart  from  approximation  theory  the  earher  results  have  obvious  analogs 
which  it  is  unnecessary  to  spell  out  in  detail.  The  new  approximation  results  corre¬ 
spond  to  lemma  2  and  theorem  5. 

Lemma  3  Assume  that  u  G  Then  we  have 

L,(l)_y(2)  <C/i^|u|2j^. 

!  W  ’ 
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Proof.  As  previously, let  Kj  denote  the  polyhedron  associated  with  the  face  i^j  and 
the  dual  edge  cr'.  By  a  Sobolev  embedding  theorem  we  have 

r  =  l,2, 


where  k]  =  ct'  and  Hence  F  :  u  is  a  bounded  linear  functional 

on  For  uniform  meshes  it  is  easy  to  verify  that  u^p  -  vanishes  for  linear 

u  so  it  follows  that 


H  2 
u)  >-u)’ 


and  the  usual  rescaling  argument  gives  that  the  right  side  is 


Then  it  follows  that 


<  Ch*  E  l-il-/ 

ixjeQ 

w'hich  gives  the  result.  ■ 

Theorem  6  Assume  that  u  €  then  we  have 

u  -  <  Ch^\u' 


w 


2,n- 


Proof.  The  proof  is  similar  to  the  proof  of  theorem  5.  ■ 


6  Numerical  test.  To  confirm  the  rate  of  convergence  given  by  the  previous 
theorem  we  computed  a  numerical  example.  The  domain  of  computation  is  the  unit 
cube.  The  domain  was  first  divided  into  equal  small  cubes  with  dimension  hxhxh. 
Then  a  dual  mesh  was  generated  (dual  cubes).  We  consider  the  following  problem 


div  u  =  0 

curl  u  =  a; 


lx=0 

1 37=1 
Uy  |j,  =  l 
Uz  \z=l 


Uy  |y=0 —  Uz  IxrrrO —  b 
sin(l)cos(t/)cos(z) 

— 2cos(x)sin(  1  )cos(2) 
cos(a:)cos(?/)sin(l) 
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where  \i  =  {Ux^Uy.  %)  and 


(cos(,T)sin(^)sin(2:)  '' 
— 2siii(a;)cos(?/)sin(2) 
sin(.T)siTi(y)cos(ir)  y 


The  exact  solution  of  this  problem  is 


u  = 


^  shi(a;)cos(y)cos(2) 
-2cos(a;)sin(?/)cos(2:) 
y  cos(.'r)cos(^)sin(z) 


Four  meshes  were  used  in  the  computation,  h  =  0.5  for  the  coarsest  mesh  and 
h  =  0.0625  for  the  finest  mesh.  The  results  are  shown  in  the  table. 


Table  1.  Errors  Between  u  and 


h 

0.5 

0.25 

0.125 

0.0625 

0.26d-01 

0.56d-02 

0.1.3d-02 

0.31d-03 

The  average  rate  of  convergence  for  this  example  is  about  which  is  slightly 
better  than  the  rate  given  by  the  theorem. 


Conclusions.  We  have  presented  an  algorithm  for  discretizing  three  dimensional 
div-curl  systems  in  three  dimensions.  We  proved  uniqueness  of  the  solution  to 
the  discretized  equations  and  an  error  estimate  showing  first  order  convergence  in 
general  and  second  order  for  regular  meshes.  A  basic  numerical  example  was  also 
shown. 
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